# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2020/11/24 10:35
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : chapter2_3.py
@Version     : Version 1.0.0
@Description : TODO
@Created By  : PyCharm
"""
import numpy as np
import matplotlib.pyplot as plt
import math
from matplotlib.font_manager import FontProperties

font_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=12)


# 求牛顿n次均差
def qiujuncha(x, f, n):
    for i in range(1, n):
        for j in range(4, i - 1, -1):
            f[j] = (f[j] - f[j - 1]) / (x[j] - x[j - i])


# 根据多项式表达式，x求y
def function(a, x, n):
    y = 0
    tmp = 1
    for i in range(n):
        y = y + a[i] * tmp
        tmp = tmp * x
    return y


# 求x的导数
def f_dao(x):
    if x == 0:
        return 0
    return 1 / (2 * math.sqrt(x))


def draw_pic_niudun(x, s):
    x1 = np.linspace(0, 64, 100)
    plt.plot(x1, lagrange(x1, x, s))
    plt.xlabel(u"x轴", fontproperties=font_set)
    plt.ylabel(u"y轴", fontproperties=font_set)
    plt.show()


# 求每一个yk/（xk-x0)(xk-x1)....(xk-xn)
def qiu_s(s, x, y):
    n = len(x)
    for k in range(n):
        tmp = 1
        for i in range(k):
            tmp = tmp * (x[k] - x[i])
        for i in range(k + 1, n):
            tmp = tmp * (x[k] - x[i])
        s.insert(k, y[k] / tmp)


# 根据拉格朗日插值函数求y
def lagrange(x1, x, s):
    n = len(x)
    result = 0
    for i in range(n):
        tmp = 1
        for j in range(0, i):
            tmp = tmp * (x1 - x[j])
        for j in range(i + 1, n):
            tmp = tmp * (x1 - x[j])
        result += tmp * s[i]
    return result


# 第一边界条件下的三次样条插值 求M
def qiu_m(h, f, o, u, d, x):
    n = len(h)
    o[0] = 1
    u[n] = 1
    d[0] = 6 * (f[0] - f_dao(x[0])) / h[0]
    d[n] = 6 * (f_dao(x[n] - f[n - 1])) / h[n - 1]
    a = []
    for i in range(1, n - 1):
        u[i] = h[i - 1] / (h[i - 1] + h[i])
    for i in range(1, n):
        o[i] = h[i] / (h[i - 1] + h[i])
    for i in range(1, n):
        d[i] = 6 * (f[i + 1] - f[i]) / (h[i - 1] + h[i])
    t = [0 for i in range(n + 1)]
    t[0] = 2
    t[1] = o[0]
    a.append(t)
    for i in range(1, n):
        t = [0 for i in range(n + 1)]
        t[i - 1] = u[i + 1]
        t[i] = 2
        t[i + 1] = o[i + 1]
        a.append(t)
    t = [0 for i in range(n + 1)]
    t[n - 1] = u[n]
    t[n] = 2
    a.append(t)
    tmp = np.linalg.solve(np.array(a), np.array(d))
    m = []
    for i in range(n + 1):
        m.append(tmp[i])
    return m


# 根据三次条插值函数求值
def yangtiao(x1, m, x, y, h, j):
    return m[j] * (x[j + 1] - x1) ** 3 / (6 * h[j]) + m[j + 1] * (x1 - x[j]) ** 3 / (6 * h[j]) + (
                y[j] - m[j] * h[j] ** 2 / 6) * (x[j + 1] - x1) / h[j] + (y[j + 1] - m[j + 1] * h[j] ** 2 / 6) * (
                       x1 - x[j]) / h[j]


def draw_yangtiao(m, x, y, h, n):
    for i in range(n - 1):
        x1 = np.linspace(x[i], x[i + 1], 10)
        plt.plot(x1, yangtiao(x1, m, x, y, h, i), color='red', label=u"三次样条插值")
    plt.xlabel(u"x轴", fontproperties=font_set)
    plt.ylabel(u"y轴", fontproperties=font_set)
    plt.title(u"三次样条插值", fontproperties=font_set)
    plt.show()


def main():
    x = [0, 1, 4, 9, 16, 25, 36, 49, 64]
    y = [0, 1, 2, 3, 4, 5, 6, 7, 8]
    n = len(x)
    s = []
    h = []
    f = y[:]
    for i in range(n - 1):
        h.append(x[i + 1] - x[i])
    u = [0 for n in range(len(h) + 1)]
    d = [0 for n in range(len(h) + 1)]
    o = [0 for n in range(len(h) + 1)]
    qiu_s(s, x, y)
    qiujuncha(x, f, 2)
    m = qiu_m(h, f, o, u, d, x)
    # print(m)
    draw_yangtiao(m, x, y, h, n)


main()
